Field of the Invention
The present invention relates to the field of imaging analysis, in particular to methods for analyzing digital images of clinical samples in a computerized method to determine boundaries of a feature of interest, such as a lesion in an MRI or CT scan. The invention also has non-medical applications.
Related Art
Presented below is background information on certain aspects of the present invention as they may relate to technical features referred to in the detailed description, but not necessarily described in detail. That is, individual compositions or methods used in the present invention may be described in greater detail in the publications and patents discussed below, which may provide further guidance to those skilled in the art for making or using certain aspects of the present invention as claimed. The discussion below should not be construed as an admission as to the relevance or the prior art effect of the patents or publications described.
In computer vision, image segmentation involves automated outlining of the boundaries of the object in the image. The goal of segmentation is to simplify and/or change the representation of an image into something that is more meaningful and easier to analyze. Image segmentation is typically used to locate objects and boundaries (lines, curves, etc.) in images. More precisely, image segmentation is the process of assigning a label to every pixel in an image such that pixels with the same label share certain characteristics.
The result of image segmentation is a set of segments that collectively cover the features of interest, or a set of contours extracted from the image (see edge/region detection). Each of the pixels in a region is similar with respect to one or more spatial characteristic or computed property, such as color, intensity, or texture. Adjacent regions are significantly similar with respect to the same characteristic(s). When applied to a stack of images, typical in medical imaging, the resulting contours after image segmentation can be used to create 3D reconstructions with the help of interpolation algorithms like marching cubes. (https (colon slash slash-en.wikipedia[dot]org/wiki/Image_segmentation #Level_set_methods)
Deformable models are popular methods for medical image segmentation and are widely used in curve evolution applications. They supply a continuous boundary of an object and are able to handle shape variations, image noise, heterogeneity, and discontinuous object boundaries. These challenging characteristics are very common in medical images. Level set models follow a non-parametric deformable model, thus are better at handling topological changes during evolution, compared to snake models. According to the mentioned above, level set can supply appropriate lesion segmentation in medical images. They can be classified into edge-based and region-based models. Edge-based models do not have a constraint of homogeneous image intensities, but are not ideal for noisy images, objects with incomplete boundaries, or for objects with low object-background contrast.
Edge-based models are not ideal for noisy images, objects with incomplete boundaries, or objects with low object-to-background contrast. Region-based models estimate spatial statistics of image regions to find the minimal energy where the model best fits the image. Chan and Vese [1] applied a region-based segmentation model with global constraint, based on the Mumford and Shah functional. The main advantage of the global constraint is its high robustness to the location of the initial contour [1,2]. However, in many cases such as with heterogeneous intensity areas, a local framework performs better than a global one [3-9]. Hybrid models are superior, defining an energy functional with local and global constraints to obtain a more accurate segmentation that is more robust to contour initialization [10]. Those methods allow curve deformation to find only significant local minima and delineate object borders despite noise, poor edge information, and heterogeneous intensity profiles [11]. Robustness to initial conditions is also an important measure of segmentation performance. Initial shape and position parameters (size, rotation, and location) need to be adequately determined; otherwise the contour may converge into a local minimum and may fail to capture the features of interest. The most common techniques for contour initialization are 1) manual selection of initial points; 2) analyzing the external force field; 3) naive geometric models such as a circle in 2-D or sphere in 3-D; 4) learned shape priors, where a statistical shape model is estimated, then the automated procedure tries to find the segmentation that best fits the shape model. However, shape prior-based methods may be quite restrictive in applications involving highly variable shapes.
Though the accuracy of contour initialization is important, the size of the local window surrounding each contour point plays a key role in the segmentation performance. The window size defines how the local scale of the statistics evaluation, and thus must be selected appropriately, even when initialization of the active contour is relatively accurate. Furthermore, well-defined local window can compensate on low-quality initial contour. Most local segmentation methods use candidates for pre-defined window sizes as input. Each candidate window size is tested over an extensive sequence of images to ascertain the best window size and this fixed window size is used for the entire database of images. However, this window size will not be optimal for all images and new images with different spatial statistics may require additional experiments to find the best window size. Thus, choosing a fixed window size by trial and error is a time consuming process. Moreover, when the images contain substantial diversity of spatial characteristics, pre-defining a single window size may result in non-optimal segmentation performance for all images. For that reason, a varied window size that is defined adaptively according to spatial information has a greater chance of providing accurate segmentation.
Some methods apply multi-scale segmentation by examining several Gaussian scales. The most accurate segmentation was obtained by using the smallest scale while with larger scale, the segmentation was more robust to contour initialization. However, for multi-scale methods a pre-specified pyramid of discrete Gaussian scales should be supplied as input. This may lead to high dependence of the segmentation accuracy on the number and the values of the discrete scales.
Recently, some works present a method to select the best scale from a range of input scales. However, choice of a single scale has some drawbacks. First, it may be sensitive to the criterion for choosing the best scale. Second, the scale choice is done by examining a specific scale each time, thus it depends on the local window only. However, in many images, global spatial characteristics can also contribute to the segmentation performance, and thus they need to be considered. Finally, these methods use the same size of the Gaussian kernel for both X and Y dimensions to define the local window (region), which may not be the best solution in cases for which image textures along the X and Y dimensions are different.
An et al. [12] implement a local framework at two different scales, which allows high sensitivity to different features. The authors show that the segmentation performance is better when using more than one single scale. Li et al. studied in depth the selection of kernel functions and its effect on the segmentation performance [5]. The authors applied their method using 3 different predefined Gaussian scales. As was expected, the most accurate segmentation was obtained by using the smallest scale while with larger scale was more robust to initialization. Yang and Boukerroui present a Gaussian scale selection based on the Intersection of Confidence Intervals rule [13]. The local scale is estimated, in the sense of minimizing the mean square error of a Local Polynomials Approximation (LPA). Pivano and Papadopoulo also applied a given pyramid of Gaussian kernels to compute local means and variances [14]. The chosen scale is the smallest one inducing an evolution speed superior to a threshold that is given by the user.
An inappropriate choice of level set parameters may lead to an inferior segmentation regardless of initialization. We propose a significant improvement of the level set segmentation method. We present an adaptive method to estimate the parameters for the level set energy function separately for each case and over iterations. When combined with estimation of an adaptive window size surrounding each contour point as suggested before, we supply a generalization of the segmentation process, applying the same model equations and deep learning architecture for any given dataset. Our method is a multi-stage process. First, we provide a novel method to estimate the parameters of the energy functional. A convolutional neural network (CNN) is used to identify the location of the zero level set contour in relation to the lesion. The output probabilities of the CNN are then used to calculate the level set parameters. Second, the adaptive window size is re-estimated by an iterative process that considers the scale of the lesion, local and global texture statistics, and minimization of the cost function over iterations. Our method requires only a single input point representing the approximate center of the detected lesion. There is no need of a more accurate initial contour as is typically supplied, automatically or manually, for local analysis. Contrary to current local active contour frameworks, our method has little to no dependence on accurate initialization and does not include any assumptions about lesion characteristics. Thus, it may perform well with highly diverse datasets that include low contrast, noisy, and heterogeneous lesions.
Few works provide an algorithm to estimate parameters for the level set energy functional. Li et al. [15] propose a fuzzy level set algorithm in which the level set parameters are initialized using a spatial fuzzy clustering approach. However, the parameters are only evaluated at the beginning of the segmentation process and remain constant throughout the whole process. In addition, the performance quality of fuzzy C-means is sensitive to noise, resulting in generally poorer segmentation. Oliviera et al. [16] present a solution for liver segmentation based on a deformable model, in which parameters are adjusted via a genetic algorithm. The genetic algorithm was used to choose the best combination of parameters from analysis of the training set, but all segmentations in the test dataset were conducted using the same parameters. Thus, this method may not be appropriate for highly diverse types of lesions. The authors in [16] also made two assumptions in their analysis: 1) the initialization is reasonably accurate; and 2) the liver is spatially homogeneous. Moreover, the authors use their method to segment the liver itself rather than a lesion dataset. The diversity of screened organs is typically much lower than the diversity that characterizes lesions. Baillard et al. [17] define the problem of parameter tuning as a classification of each point along the contour. That is, if a point belongs to the object, then the surface should locally extend, and if not, the surface should locally contract.
This classification is performed by maximizing the posterior segmentation probability [17]. However, the authors only consider the direction of the curve evolution and not its magnitude, which is critical especially inheterogeneous regions, wherein convergence into local minima should be prevented. Thus, both [16] and [17] are likely to have limited performance for highly diverse datasets, given the limited amount of information that is incorporated. We evaluated our proposed adaptive framework with two different energy models: the piecewise constant model and the mean separation model. Details of each model follow.
A. Piecewise Constant Model (PC)
Piecewise constant models are region-based models that assume intensity homogeneity. Chan and Vese [1] present a global framework of this PC model, which assumes intensity homogeneity. Let FPC(Mu, Mv, ϕ) be a function which models the object and its background as uniform intensities represented by their means Mu and Mv:
                                          F            PC                    ⁡                      (                                          M                u                            ,                              M                v                            ,              ϕ                        )                          =                              μ            ⁢                                          ∫                Ω                            ⁢                                                δϕ                  ⁡                                      (                                          x                      ,                      y                                        )                                                  ⁢                                                                        ∇                                          ϕ                      ⁡                                              (                                                  x                          ,                          y                                                )                                                                                                              ⁢                dxdy                                              +                                    λ              1                        ⁢                                          ∫                Ω                            ⁢                                                                                                                                      I                        ⁡                                                  (                                                      x                            ,                            y                                                    )                                                                    -                                              M                        u                                                                                                  2                                ⁢                H                ⁢                                                                  ⁢                                  ϕ                  ⁡                                      (                                          x                      ,                      y                                        )                                                  ⁢                dxdy                                              +                                    λ              2                        ⁢                                          ∫                Ω                            ⁢                                                                                                                                      I                        ⁡                                                  (                                                      x                            ,                            y                                                    )                                                                    -                                              M                        v                                                                                                  2                                ⁢                1                                              -                      H            ⁢                                                  ⁢                          ϕ              ⁡                              (                                  x                  ,                  y                                )                                      ⁢            dxdy                                              (        1        )            where μ affects the smoothness of the curve, and (x,y) is an examined location within a given image I. Ω is a bounded subset in R2 and ϕ(x,y) is a signed distance map. ∇ is the first variation of the energy with respect to the distance map ϕ(x,y). Let C be a parameterized curve, closed contour in Ω represented by the zero level set (ZLS). C={(x,y)|ϕ(x,y)=0}. The interior region of C is defined by the smoothed Heaviside function Hϕ(x,y).
                              H          ⁢                                          ⁢                      ϕ            ⁡                          (                              x                ,                y                            )                                      =                  {                                                                      1                  ,                                                                                                  ϕ                    ⁡                                          (                                              x                        ,                        y                                            )                                                        >                  ɛ                                                                                                      0                  ,                                                                                                  ϕ                    ⁡                                          (                                              x                        ,                        y                                            )                                                        <                                      -                    ɛ                                                                                                                                            1                    2                                    ⁢                                      {                                          1                      +                                                                        ϕ                          ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                          ɛ                                            +                                                                        1                          π                                                ⁢                                                  sin                          ⁡                                                      (                                                                                          πϕ                                ⁡                                                                  (                                                                      x                                    ,                                    y                                                                    )                                                                                            ɛ                                                        )                                                                                                                }                                                                                                                                                              ϕ                      ⁡                                              (                                                  x                          ,                          y                                                )                                                                                                  <                  ɛ                                                                                        (        2        )            The external region of C is defined by (1-Hϕ(x,y)). Representation of the narrow band around the ZLS contour C can be done by using the derivative of Hϕ(x,y) and obtaining a smooth approximate version of Dirac delta δϕ(x,y).
By applying a local version of the PC model, Mu and Mv can be replaced by their local equivalent terms, mu and mv. In that case, they will represent the local means of a region surrounding each contour point [3].
B. Mean Separation (MS) Model
The Mean Separation model was first proposed by Yezzi et al. [2]. It assumes that object and its background should have maximally separate mean intensities:
                                          F            MS                    ⁡                      (                                          m                u                            ,                              m                v                            ,              ϕ                        )                          =                              μ            ⁢                                          ∫                                  Ω                  n                                            ⁢                                                δϕ                  ⁡                                      (                                          x                      ,                      y                                        )                                                  ⁢                                                                        ∇                                          ϕ                      ⁡                                              (                                                  x                          ,                          y                                                )                                                                                                              ⁢                dxdy                                              +                                    λ              1                        ⁢                                          ∫                                  Ω                  n                                            ⁢                                                                                                                                                                I                          ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                          -                                                  m                          u                                                                                            A                        u                                                                                                  2                                ⁢                H                ⁢                                                                  ⁢                                  ϕ                  ⁡                                      (                                          x                      ,                      y                                        )                                                  ⁢                dxdy                                              +                                    λ              2                        ⁢                                          ∫                                  Ω                  n                                            ⁢                                                                                                                                                                I                          ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                          -                                                  m                          v                                                                                            A                        v                                                                                                  2                                ⁢                                  (                                      1                    -                                          H                      ⁢                                                                                          ⁢                                              ϕ                        ⁡                                                  (                                                      x                            ,                            y                                                    )                                                                                                      )                                ⁢                dxdy                                                                        (        3        )            
Ωn∈R2 is the local version of Ω which represents the narrow-band points, each point with its local surrounding only. Au and Av are the areas of the local interior and exterior regions surrounding a contour point, respectively:
                                          A            u                    =                                    ∫                              Ω                n                                      ⁢                          H              ⁢                                                          ⁢                              ϕ                ⁡                                  (                                      x                    ,                    y                                    )                                            ⁢              dxdy                                      ,                              A            v                    =                                    ∫                              Ω                n                                      ⁢                                          (                                  1                  -                                      H                    ⁢                                                                                  ⁢                                          ϕ                      ⁡                                              (                                                  x                          ,                          y                                                )                                                                                            )                            ⁢              dxdy                                                          (        4        )            When mu and mv are the most different from each other, a minimization of the MS energy can be obtained. In some cases, the MS model may able to supply better results than the PC model due to preference of maximal contrast between the interior and the exterior regions without any restrictions on regions homogeneity. This allows to find image edges very well without considering if the interior or exterior regions are uniform or not.C. Histogram Separation (HS) Model
Consider Pu(b) and Pv(b) to be two intensity histograms computed from the local regions, interior and exterior respectively, surrounding each ZLS contour point. N is the number of histogram bins. The Bhattacharyya coefficient B is chosen to evaluate the similarity of the two histograms [18]:
                    B        =                              ∫                          b              =              1                        N                    ⁢                                                                                          P                    u                                    ⁡                                      (                    b                    )                                                  ⁢                                                      P                    v                                    ⁡                                      (                    b                    )                                                                        ⁢            db                                              (        5        )            By separating the intensity histograms, the interior and exterior regions are allowed to be heterogeneous as long as their whole intensity profiles are different. By using the Bhattacharyya metric to quantify the separation of intensity histograms, we are able to segment objects that have local non-uniform intensities. There is no preliminary assumption regarding the gray level distribution of the object. However, the intensity profile of the entire object and the entire background must still be separable.
                                          F            HS                    ⁡                      (                          x              ,              y              ,              ϕ                        )                          =                                            ∫                              Ω                n                                      ⁢                                          B                ⁡                                  (                                                            1                                              A                        u                                                              -                                          1                                              A                        v                                                                              )                                            ⁢              dxdy                                +                                    λ              1                        ⁢                                          ∫                                  Ω                  n                                            ⁢                                                (                                                            1                                              A                        u                                                              ×                                                                                                                        P                            u                                                    ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                                                                                P                            v                                                    ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                                                                ⁢                    H                    ⁢                                                                                  ⁢                                          ϕ                      ⁡                                              (                                                  x                          ,                          y                                                )                                                                              )                                ⁢                dxdy                                              -                                    λ              2                        ⁢                                          ∫                                  Ω                  n                                            ⁢                                                (                                                            1                                              A                        v                                                              ×                                                                                                                        P                            v                                                    ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                                                                                P                            u                                                    ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                                                                ⁢                                          (                                              1                        -                                                  H                          ⁢                                                                                                          ⁢                                                      ϕ                            ⁡                                                          (                                                              x                                ,                                y                                                            )                                                                                                                          )                                                        )                                ⁢                dxdy                                                                        (        6        )            